function [p,e,t] = new2d_grid(N,limits)
% [p,e,t] = new2d_grid(N,limits)
%
% Build a 2D structured grid using N(1) subdivisions in x and N(2)
% subdivisions in y. The domain is the square
% [limits(1),limits(2)]x[limits(3),limits(4)].
% If all the N(:) are even, it is guaranteed that each element has at
% least one internal node.

 nv = (N(1)+1)*(N(2)+1);
 p = zeros(2,nv);
 for i=1:N(1)+1
   for j=1:N(2)+1
     ij = (i-1)*(N(2)+1) + j;
     p(:,ij) = [i,j]' - 1;
   end
 end
 for i=1:2
   j = (i-1)*2 + 1;
   k = j+1;
   p(i,:) = ((limits(k)-limits(j))/N(i))*p(i,:) + limits(j);
 end

 ne = 2*N(1)*N(2);
 t = zeros(3,ne);
 e = [];
 for i=1:N(1)
   if(2*floor(i/2)==i)
     oj = 'south';
   else
     oj = 'north';
   end
   for j=1:N(2)
     ij = (i-1)*N(2)*2 + (j-1)*2;
     iv(1) = (i-1)*(N(2)+1) + j;
     iv(2) = (i-1)*(N(2)+1) + j+1;
     iv(3) = (i-1+1)*(N(2)+1) + j;
     iv(4) = (i-1+1)*(N(2)+1) + j+1;
     [tt,ex,ey] = new2d_grid_singlesquare(iv,oj);
     t(:,ij+1:ij+2) = tt;
     if(i==1)
       e = [e,ex(:,:,1)];
     endif
     if(i==N(1))
       e = [e,ex(:,:,2)];
     endif
     if(j==1)
       e = [e,ey(:,:,1)];
     endif
     if(j==N(2))
       e = [e,ey(:,:,2)];
     endif
     if(strcmp(oj,'north')==1)
       oj = 'south';
     else
       oj = 'north';
     endif
   end
 end

return

